gpu = gpuDevice();
sizes = power(2, 20:2:23);
repeats=5;
N = sqrt(sizes);
mmTimesHost = inf(size(sizes));
mmTimesGPU = inf(size(sizes));
for ii=1:numel(sizes)
    A = rand( N(ii), N(ii) );
    B = rand( N(ii), N(ii) );
    % First do it on the host
    for rr=1:repeats
        timer = tic();
        C = A*B; %#ok<NASGU>
        mmTimesHost(ii) = min( mmTimesHost(ii), toc(timer));
    end
    % Now on the GPU
    A = gpuArray(A);
    B = gpuArray(B);
    for rr=1:repeats
        timer = tic();
        C = A*B; %#ok<NASGU>
        wait(gpu);
        mmTimesGPU(ii) = min( mmTimesGPU(ii), toc(timer));
    end
end
mmGFlopsHost = (2*N.^3 - N.^2)./mmTimesHost/1e9;
[maxGFlopsHost,maxGFlopsHostIdx] = max(mmGFlopsHost);
mmGFlopsGPU = (2*N.^3 - N.^2)./mmTimesGPU/1e9;
[maxGFlopsGPU,maxGFlopsGPUIdx] = max(mmGFlopsGPU);
fprintf('Peak calculation rate: %1.1f GFLOPS (host), %1.1f GFLOPS (GPU)\n',  maxGFlopsHost, maxGFlopsGPU)